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Abstract 

A fluid between two spheres, concentric or not, at different temperatures will flow in the presence 
of a constant gravitational force. Although there is no possible hydrostatic state, energy transport is 
dominated by diffusion if temperature difference between the spheres is small enough. In this conductive 
regime the average Nusselt number remains approximately constant for all Rayleigh numbers below some 
critical value. Above the critical Rayleigh number, plumes appear and thermal convection takes place. 
We study this phenomenon, in particular the case where the inner sphere is displaced from the centre, 
using a two-component thermal lattice Boltzmann method to characterize the convective instability, the 
evolution of the flow patterns and the dependence of the Nusselt number on the Rayleigh number beyond 
the transition. 


1 Introduction 

Natural convective flow between a sphere at constant temperature and its spherical enclosure at a different 
temperature is an idealisation of many problems of practical interest. Convection patterns between concentric 
spherical shells were first observed by Bishop et al. in 1966 for a range of aspect ratios [1]. Experiments 
were performed with the inner sphere hotter than the outer sphere. A wide variety of steady and non steady 
patterns were observed in subsequent experiments using air, water and silicone oil; including the case of 
vertically eccentric spherical shells [2]. For concentric spheres, analytic solutions for the conductive regime 
were reported by Make & Hardee [3]. The dependence of the Nusselt number Nu on the Rayleigh number 
Ra on the transition from the conductive to the convective steady solutions was obtained by Teertstra et 
al. [4]. The stability analysis of the convective regime has deserved wide attention, see for example [5]. 
Numerical simulations of the concentric configuration can be found for a variety of numerical methods, in 
particular for the study of the transition to turbulence [SI El El m HQ]. 

Little is known of the inverted configuration. The case of an outer shell hotter than an inner concentric 
shell was experimentally and numerically studied by Futterer et al. using silicone oil m- They found 
unsteady periodic flows consisting of cold blobs dripping from the inner sphere at a frequency that became 
irregular through a period-doubling process as Ra increases. 

While concentric configurations have received much attention, few studies can be found about the eccen¬ 
tric ones. In the present work, we report a numerical exploration of the natural convection between eccentric 
spherical shells at fixed temperatures (see figure [^. Simulations were performed using a three dimensional, 
two-component lattice Boltzmann method (LBM) that approximates solutions to the Oberbeck-Boussinesq 
equations of thermal convection m- This method has been used successfully for the simulation of natural 
convection phenomena HSllIl]. 

In the next section we present the numerical method and its validation. We compare predicted results 
with concentric and vertically eccentric numerical and experimental observations found in the literature. We 
present numerical solutions for different eccentric configurations in the third section. The transition from 
conductive to convective heat transport is characterised, and a series of steady and unsteady patterns and 
behaviour are presented in this section. In the last section we summarise our observations and motivate 
future work. 

^Author to whom correspondence should be addressed: cmi.ciencias@ciencias.unam.mx 
Whysics Department, School of Science, Universidad Nacional Autonoma de Mexico 
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2 The thermal LBM 


To simulate natural convection phenomena we used a D3Q19 two-component lattice Boltzmann equation. 
The method is a direct three dimensional extension of that proposed by Innamuro in 2002 m- Space is 
discretized using a cubic lattice where a density fk and a temperature gf. distribution functions are computed. 

The distribution functions are then used to compute the fluid velocity u and temperature T at the lattice 

nodes. Lattice spacing as well as time steps can be conveniently set to unity. At every node r in the lattice, 
the distribution functions evolve in time according to 

fk{r + ek,t + l) = fk{r,t) - \[fk{r,t) - fl'^{r,t)]+Gk, (1) 

gk{r + ek,t + l) = gk{r,t)-^^[gk{r,t)-gl'^{r,t)]. ( 2 ) 

The coefficients r and Tg represent relaxation times and are related to the fluid kinematic viscosity 
u = {r — l/2)/3 and thermal diffusivity a = (rg — l/2)/3. The local equilibrium distribution functions 
and are given by 
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The equilibrium distributions depend on the macroscopic fields n, T and p, the mass density, and must 
be computed every time step through 
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The constants Wk in the equilibrium functions definitions and 0 take the values wq = 1/3, = 1/18 

for k = 1, ...,6 and = 1/36 for k = 7,..., 18. The last term in § is related to the body force and gives the 
buoyancy term in the Boussinesq equations. It is defined as Gk = —3/3wk {T{r,t) — Tq) where (3 is the 
coefficient of thermal expansion of the fluid, g the acceleration due to gravity and Tq a reference temperature 
taken as the average of the temperatures of the inner and the outer shells. 

The set of microscopic velocities {ck : /c = 0,..., 18} is given by 


eo = (0,0,0), Cl = -64 = (1,0,0), 62 = -65 = (0, 1,0), 63 = -66 = (0,0, 1), 

67 = -610 = (1,1,0), eg = -611 = (1,0, 1), 69 = -612 = (0, 1, 1), 

613 = —6i6 = ( —1 , 1, 0), 614 = —617 = ( —1, 0, 1 ), 615 = —6i8 = (0, —1 , 1). (8) 

Equations 0 and Q, with the choice microscopic velocities, provide an algorithm for updating all the 
distribution functions fk and gk at a given node in the lattice, as long as its 18 nearest neighbours in 
the lattice are inside the fluid domain. For nodes close enough to a solid wall, the distribution functions 
coming from neighbouring nodes outside the fluid domain must be provided as a boundary condition for the 
method. We choose to adopt the set of boundary conditions proposed by mi for curved solid walls at fixed 
temperatures. In this boundary conditions, distribution functions that should provide the node outside the 
fluid domain are extrapolated considering that T and u are prescribed at the point in the boundary that 
intersects the line joining the nodes in question. 

The method was implemented to run in parallel in Graphic Processing Units (GPU) because of the large 
number of nodes needed in the simulation of three dimensional flows. Typically, around 10^ nodes where 
used to obtain unsteady flows and simulations took a few days running on a single GPU. 
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Figure 1 : (a) Schematic representation of the problem, (b) Steady numerical solution for 77 = 0.5, ^ = tt, 

e = 0.6 and Ra = 1.875 x 10^. The numerical solutions is symmetric with respect to the vertical axis and 
the figure shows isothermal curves in the left half of a plane of symmetry and stream lines in the right half. 


3 Problem statement and validation 


A fluid of initial density po fills the gap between two eccentric spherical shells at fixed positions. The inner 
sphere of radius is at temperature and the external sphere of radius rg is at a higher temperature Tg. 
The position of the centre of the inner sphere is given by the polar angle 0 and the distance to the external 
sphere centre c, as shown in figure [^a). The vertical plain containing the centres of both spheres is a plane 
of symmetry of the fluid domain. 

The thermal LBM presented in the previous section approximates solutions to the Boussinesq equations 
describing the convective flow of the stated problem [ 12 ]. Scaling lengths with L = rg — temperatures 
with AT = T^ — Tg, velocities with a/L and pressure with poo^/L‘^, the convective flow is defined by the 
non-dimensional eccentricity e, the polar angle 6 >, the aspect ratio 77 , the Rayleigh number Ra and the Prandtl 
number Pr given by 



ua 



( 9 ) 


For all the flows studied Pr = 0.7, corresponding to air, and T^ > Tg which gives Ra > 0. 

To validate the code we compared with a variety of results in the literature. Numerical results reproduced 
the transition from conductive to convective energy transport observed experimentally for concentric spheres 
by Teertstra et al. [4] through the measurement of the average Nusselt number on the internal shell, defined 
here as 


Nu = 
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where n is the normal direction to the spheres. 

The values of Nu obtained for steady axisymmetric flows are compared to those of Teertstra et al in 
figure]^ for different values of 77 , showing a good agreement with experiments. The insets in figure [^a) show 
isothermal curves (left half) and stream lines (right half) obtained with the code on a plane of symmetry. 
Each inset shows the obtained flow patterns and plume formation through three characteristic regions, 
roughly divided by the vertical dotted lines. These regions correspond to the conductive and convective 
regimes, and a transitional region in between them. 

A linear fit of the numerical results for the steady convective region in the range 0.2 < 77 < 0.7 suggests 
a power law relation given by 

Nu = 0.3Q{l-r])Ra°-^^'^+°-'^^. ( 11 ) 


This relation is consistent with correlations obtained by Raithby & Hollands m, Scanlan et al. m and 
Feldman & Colonius | 6 ]. 

In figure [^b) isothermal curves and stream lines are shown for a case of vertically eccentric spheres. The 
flow pattern is in good agreement to that observed experimentally by Powe et al. [ 2 ] for the set of parameters 
in the region they named “Steady Multiple Interior Cells”. 
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Figure 2 : Comparison of the steady state Nu numbers obtained with LBM and experiments performed by 
Teertstra et al. [4] for concentric shell and different aspect ratios. 77 = 1/3 in (a), 77 = 5/24 in (b), 77 = 1/2 
in (c) and 77 = 2/3 in (d). Insets in (a) show typical isothermal lines in the left half, and steam lines in the 
right half, found in the three regions divided by the dotted vertical lines. 


Convection between concentric spheres for 77 = 0.714 served to validate the code with unsteady three 
dimensional flows. Figure shows numerical solutions for Ra = 4.03 x 10^, just above the critical value 
{Ra = 3.96 X 10^) predicted by the linear stability analysis presented by Travnikov et al. [5]. The non- 
dimensional temperature 0 = {T — Ti)/AT and azimuthal velocity component on the midrange spherical 
surface r = {ve T Tj)/2 show an azimuthal mode m = 10, which is also the critical wave number predicted 
by linear stability. This unsteady periodical flow was obtained by letting the solution reach a steady state 
at Ra = 3.8 x 10^ before increasing its value to Ra = 4.03 x 10^. The pattern of u^f) is very similar to that 
obtained by Scurtu et al. [7] for travelling and pulsating waves of mode ttx = 10 using spectral methods, and 
those obtained by Feldman & Colonius [6] for mode 772 = 12 at higher Ra values using openFoam algorithms. 
When Ra was increased to 5 x 10^, the solution showed an ttx = 11 mode over many oscillations before 
changing to an ttx = 12 mode, which is consistent with results in mi]. The non-dimensional velocity 
shown in figure]^ rescaled taking a = 2 x 10~^m‘^/s and L = O.Ittx, predicts azimuthal velocities of the 
order of 0.1m/s. 

In general, steady and time averaged Nu numbers predicted by the code are in good agreement with 
the corresponding values previously reported for steady and periodic flows. This is summarised in Table 
where numbers with an asterisk superscript denote a time average Nu of a periodic solution. Notice that 
a steady state was obtained for 77 = 0.667 and Ra = 10^, this is in agreement with experiments in [4] (see 
figure l^c)) . 

Typical simulations required 4 x 10^ time steps and the three dimensional fluid domain was represented 
by as much as 2 x 10^ grid points for unsteady solutions. Using a Tesla C2075 parallel processor, runs took 
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Figure 3: Numerical solutions of the non-dimensional temperature and azimuthal velocity component for 
T] = 0.714, e = 0 and Ra = 4.03 x 10^ on the midrange spherical surface r = (re + ri)/2 showing an azimuthal 
mode m = 10. In (a) and (c) 0 is shown from the top and lateral views respectively. In (b) and (d) U(f) is 
shown from the top and lateral views respectively. 


between three and six days. 


4 Eccentric configurations 


For the case of concentric shells, the average Nu number defined in ( |10| ) represents the average heat flux 
on either the external or the internal sphere, divided by kATr]/L or kAT/{Lr]) respectively. The latter 
expressions are the average heat flux for the convection problem in concentric spheres in the external and 
internal spheres respectively. Defined in this way Nu tends to one in both spheres for small values of Ra as 
heat transport is dominated by conduction (see figure]^. For eccentric spheres we defined the average Nu 
number accordingly as 


Nu = 


L 

47rr^ATA/'?i* 



( 12 ) 


Table 1: Steady and time averaged Nu numbers 


T] 

Ra 

Present 

Ref. lU 

Ref. [g 

Ref. E] 

Ref. llflj 

Eq. ( 

ET 

0.5 

10^ 

1.0077 

1.0217 

1.001 

1.000 




0.5 

10® 

1.0865 

1.104 

1.0990 

1.1310 

1.1021 

1.0478 

0.5 

10® 

1.8977 

1.9665 

1.9730 

1.9495 

1.9110 

1.88848 

0.5 

10® 

3.3077 

3.4012 

3.4890 

3.4648 

3.355 

3.3906 

0.667 

10® 

1.0415 

1.04825 

1.001 

1.00115 




0.667 

10® 

1.7295 

1.793* 

1.073 

1.07138 


1.7886 

0.833 

10® 

1.0118 

1.011 

1.0 

1.0018 




0.833 

10® 

1.5993* 

1.6523* 

1.001 

1.0028 
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Figure 4: Numerical prediction of the steady Nu number showing the convective to conductive transition 
for eccentric configurations for 77 = 0.5 and e = 0 . 6 . Solutions at polar angles 0 = jj 
are compared with the concentric solution represented by the continuous line. The insets show the typical 
isothermal curves (top) and stream lines (bottom) for ^ = 37 r /4 on the plane of symmetry in the conductive 
and convective regions, and a transitional region in between them marked by the vertical dotted lines. 


where r takes the values Ve and and Nu^ is the average Nusselt number for the heat conduction problem 
between eccentric spheres obtained by Alassar m- Through bispherical coordinates, is expressed as 
an infinite series 


Nu* = V _ - _ 

^2 at ^ _ e(2"+i)Ci ’ 

n=0 


(13) 


where r takes the values re and r^, = smh ^(a/re), Ci = sink ^(a/r^) and 

a = v/(c + + re)(c + r^ - re)(c - u re)(c - u - re)/2c 


determines the foci of the bispherical coordinates. 

We first searched for steady solutions for fixed values of r] and e and different angles 0. In figure 
appears the behaviour of the steady Nu number as Ra increases for 77 = 0.5, e = 0.6, and the set of angles 
0 = 0 , 7 r/ 4 , 7 r/ 2 , 37 r/ 4 , 7 r. Numerical solutions show a similar behaviour of that observed in concentric spheres, 
given by the continuous line in the figure. There is a transition from a conductive to a convective regime at 
roughly the same values of the Ra number as the concentric solutions. At angles 6 > = 37 r /4 and tt the steady 
convective regime seems to be reached at lower values of Ra than the rest of cases shown in the figure. Also, 
the slopes of the convective regime are slightly lower for the eccentric configurations. A linear fit suggest 
exponents of values around 0.245 for 0 > 7 r/ 2 , 0.225 for 0 < 7 r /2 and its lowest value of 0.2 for 0 = 7 r/ 2 . 
Insets show the evolution of isothermal curves and steam lines. On the other hand, eccentric simulations 
become unstable at lower values of the Ra number when compared to concentric ones. 

Beyond the steady convective region a series of periodic solutions where found. Starting from an oscillating 
Nu number at the inner and outer spheres with a clear frequency, as the Ra number increased unsteady 
flow behaviour became irregular through what seamed as a period doubling process. The behaviour Nu 
along non-dimensional time t is shown in figure for the case of 77 = 0.5, e = 0.8 and 0 = 7 r/ 4 . Taking 
a = 2 X and L = O.Itti, the times on figure]^ must be rescaled with a time of 500s and so the 

figure would show a lapse of 10 s. 
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Figure 5: Time evolution of the average Nu number calculated on the external sphere for the case of = 0.5, 
e = 0.8 and ^ = 7r/4 and increasing values of the Ra number. The different curves show a periodic behaviour 
of increasing complexity as the Ra number grows. Details of the flows at the points indicated with the 
symbols o, o, A and V are shown in figures and 


Figure [^shows in columns some details of the unsteady convecive flow for Ra = 6.25 x 10^ at three times, 
indicated in figure when the Nu number reaches an absolute maximum value, the following minimum and 
the next local maximum respectively (the symbols o, o and A in figure [^. The top row shows the isothermal 
surface 0 = 0.7, the row in the middle shows the temperature field on the plane of symmetry of the flow 
domain and the bottom row shows a set of stream lines superposed to the surface of 0 = 0.7. 

Observations suggest that isothermal surfaces suffer an instability in the direction transversal to the 
plane of symmetry, generating a sort of “tail” at the top of the surface that seem to flap vertically. This 
oscillation is accompanied by currents coming from the “front” of the inner sphere when the Nu number 
reaches its peaks (see figures [^g) and|^i)). Lower values of Nu appear with smoother convection cells and 
flow patterns within the flow (see figure [^h)). At higher Ra numbers, the currents coming from the front 
start to oscillate horizontally showing a new instability of the isothermal surfaces. This can be seen in figure 
corresponding to the point marked with the symbol V in the curve of Ra = 6.8 x 10^ in figure 

5 Conclusions 

It is remarkable that such a simple numerical model as LBM can reproduce steady and unsteady three 
dimensional natural convection flows. Additionally, the simplicity of the algorithm allows for an efficient 
implementation on massively parallel architectures. This was the motivation for exploring natural convec¬ 
tion in eccentric annuli configurations with this method. Although there is much work done in concentric 
configurations, eccentric ones have received little attention and our observations contribute to the complete 
picture of natural convection in spherical annuli. 

When validated with numerical and experimental results of the concentric configuration, results suggested 
a power law correlation of the average Nu number on the spherical boundaries with the Ra number and 
the aspect ratio, see equation close to what has been proposed in previous work (see [m m E]). 

When simulating eccentric configurations, result suggested a transition from a conductive regime to a steady 
convective regime similar to that found in the concentric case. Convective regimes are reached at lower Ra 
numbers and Nu grows at a smaller rate in eccentric configurations. The lack of symmetry should produce 
more mixing in and lower the heat flux observed in the convective regime. 
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Figure 6: Results for the case of r] = 0.5, e = 0.8, 0 = 7r/4 and Ra = 6.8 x 10^ at the times indicated by the 
points o, o and A in figure The row at the top show surface 0 = 0.7 located inside the enclosing sphere. 
The row in the middle shows a slice of the temperature field at the plane of symmetry of the configuration. 
The bottom row shows stream lines seeded on the white line coloured with the magnitude of the velocity. 


Eccentric convective states also become unstable at lower values of Ra. As Ra grows, flow structure 
and the Nu number present periodic fluctuations in time. Nu oscillations are related to the deformation 
of isothermal surfaces. These show a modulation transversal to the vertical plane of symmetry and, as Ra 
grows, surfaces start to oscillate in the horizontal direction as well. Similar behaviour of the periodicity of 
the flow and heat flux was found in natural convection in a concentric annuli for an inner sphere colder than 
the outer one m- 

Results here presented are encouraging and suggest that many different examples of three dimensional 
natural convection phenomena between a body and its enclosure can be studied using LBM. Specially given 
the simplicity to impose boundary conditions on curved surfaces in LBM, and the efficiency that exhibits 
under fine grain parallelism. 
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Figure 7: Isothermal surface 0 = 0.7 superimposed to a set of stream lines seeded on the white line at the 
time indicated by V in figure Stream lines are coloured with the magnitude of the velocity. Parameter 
values are 77 = 0.5, e = 0.8, = 7 r /4 and Ra = 6.8 x 10^. 
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